Simple Lattice-Models of Ion Conduction: 
Counter Ion Model vs. Random Energy Model 



(N 
O 

o 

(N 

C 
3 



(N 



i 

c 

o 
o 



(N 

o 

O 
(N 
O 

03 



c 

o 
o 



X 
S3 



J. Rcinisch and A. Heuer 

Westfalische Wilhelms-Universitat Miinster 
Institut fiir Physikalische Chemie and SFB 458 
Schlossplatz 4/% 48149 Miinster, Germany 
(Dated: February 1, 2008) 

The role of Coulomb interaction between the mobile particles in ionic conductors is still under 
debate. To clarify this aspect we perform Monte Carlo simulations on two simple lattice models 
(Counter Ion Model and Random Energy Model) which contain Coulomb interaction between the 
positively charged mobile particles, moving on a static disordered energy landscape. We find that 
the nature of static disorder plays an important role if one wishes to explore the impact of Coulomb 
interaction on the microscopic dynamics. This Coulomb type interaction impedes the dynamics 
in the Random Energy Model, but enhances dynamics in the Counter Ion Model in the relevant 
parameter range. 

PACS numbers: 66.30.Dn 



I. INTRODUCTION 

Many disordered insulating materials show univer- 
sal behavior in their ionic DC- and AC-conductivity. 
Two prominent examples are the Arrhenius tempera- 
ture dependence of the DC-conductivity and the fact 
that AC-conductivity data for different temperatures can 
be scaled onto a single master curve. This curve is 
similar for most materials. This observed universality 
jl], 1^] has stimulated researches to find a common mecha- 
nism for ion conductivity in these materials. A different 
type of system are disordered Anderson insulators, dis- 
playing electron transport. They are denoted Coulomb 
glasses. In contrast to ion conductors non-localized dy- 
namical processes play an important role in these sys- 
tems. The non-localized dynamics give rise to a very dif- 
ferent temperature dependence of the DC-conductivity at 
low temperatures, which turns out to be proportional to 
expi-iTM/T) 1 / 4 } (Mott's law §) or to exp^To/I 1 ) 1 / 2 ] 
(Efros-Shklovskii law (§). 

Similar computational approaches have been chosen 
to investigate both types of problems and various the- 
oretical models have been developed in this context 
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§, |, @, §|, 0. The focus of this work lies on the ef- 
fect of Coulomb interaction on the microscopic dynamics 
in ion conductors. We chose two similar lattice models 
to investigate the relevant microscopic dynamics: The 
Counter Ion Model (CIM) Q, || and the Random En- 
ergy Model with cation-cation interaction (REM) [ p"3| |, 
which has some features in common with the Coulomb 
glass model used for the second group of disordered solids 
[g, [LCI]]. The CIM and the REM are designed to reflect 
important aspects of vitreous ion conductors, which are 
a high degree of disorder and mobile charged ions in a 
fixed glass network. The models differ in the way how the 
time independent (static) disorder is realized. In this pa- 
per we analyze the effects of Coulomb interaction among 
mobile ions and show that the nature of the disorder has 
great influence on the dynamics. In Sect. || we present 
the mode ls as well as computational details, while in the 
Sect. Ill the results are presented and discussed. 



Figure 1: One dimensional illustration of the static potential 
landscape in the random energy model. The arrows indicate 
possible jumps. 



II. MODELS AND COMPUTATIONAL DETAILS 

A. Similarities of both models 

Both models base on a single type of mobile particles 
restricted to discrete sites in a simple cubic lattice with 
I 3 sites of distance a. The distance for a single jump is 
limited to a, and occupied positions arc forbidden. All 
mobile particles have the same positive charge. In con- 
trast to other works on the CIM jll], |l2| the strength 
of the Coulomb interaction between the mobile particles 
can be varied independently from the strength of static 
disorder by appropriate selection of parameters, as it is 
also possible for the REM; see below for precise defini- 
tions. This variation is essential to elucidate the effect of 
Coulomb interaction among the mobile particles on their 
dynamics. There are two energy contributions to the 
total Hamiltonian: Mobile particles moving on a time 
dependent potential surface (generated by the particles 
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Figure 2: One dimensional illustration of the static potential 
landscape in the counter ion model. The arrows indicate the 
possible jumps. 
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themselves) and on a time independent potential surface. 
We call the first dynamic and the latter static. The dy- 
namic part has the same form in both models and the 
Hamiltonian for the cation-cation interaction is, 
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(1) 



A configuration with sites i is described by the occupa- 
tion numbers n.; = 1 for occupied sites and = for 
empty sites. The omitted factor 47reo is taken in account 
by use of appropriate units. The mean nearest-neighbor 
interaction in a system with randomly distributed cations 
can be written as 



Vc = -■ 



(2) 



r s is the mean distance of nearest neighbors in such a 
system with cation-concentration c = N/l 3 , 
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Furthermore we introduce the dimensionless parameter 
r ca t via 



Vc 



k B T k B Tr s 



(4) 



In the literature this parameter is already established for 
the REM Jl3| and we use it for the CIM to ensure compa- 
rability. The Hamiltonian for the cation interaction can 
thus be rewritten as, 



k B T 



V 



(5) 



In what follows V c is regarded as a constant for a given 
concentration. Here r cat oc e 2 is a measure for the inter- 
action strength between the cations relative to k B T. 



B. Static Disorder in the Random Energy Model 

The REM features a straightforward approach of 
the principle of representing complexity by randomness. 



Figure 3: Typical errors and finite size effects for both models. 
The largest errors are of the size of the symbols and occur at 
large t. 



Each site i gets a random energy taken from a Gaussian 
distribution with standard deviation a e and mean value 
0. The Hamiltonian for the static disorder becomes, 



Hstat — 



E? 



(0) 



The model is fully determined by two dimensionless pa- 
rameters: a r is connected to the standard deviation of 
the random site energies by, 



(7) 
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In changing oy one can change the strength of static dis- 
order relative to temperature in the system. As before, 
the second parameter T ca t determines the cation-cation 
interaction strength. The total Hamiltonian reads, 
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In Fig. the potential landscape of the REM is illus- 
trated. 



C. Static Disorder in the Counter Ion Model 

In the CIM static disorder is generated by randomly 
placing negatively charged particles at centers of cubic 
lattice cells. As for cations these anions cannot occupy 
the same site. The resulting potential landscape (see Fig. 
B) is different from that seen in Fig. ^ As for the REM 
we used T cat to control the interaction strength among 
cations and defined a parameter r s t at adequate to control 
the interaction strength between cations and anions. It 
is defined in analogy to r cat as, 



F stat — 



k C a ' V c 

k B T 



(9) 



Figure 4: Left: Dependency of the low frequency diffusion 
constant Ddc on ay in the REM. The lines correspond to 
exponential fits, the data for o> = are omitted. The dot- 
ted line describes the temperature dependency of Ddc for a 
constant quotient j? 21 - = -pf- = 0.05 . The errors are smaller 
than the symbol size. 

Right: Dependency of the slope on F ca t- The data match the 
theoretical limit at T ca t = 0. 
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Figure 5: Dependency of the high frequency diffusion constant 
Dac on ay in the REM. The lines are exponential fits and 
serve as guides to the eye. The errors are within the symbol 
size. 



III. RESULTS 



The factor k ca has been introduced to vary the strength of 
the cation-anion interaction separately from the cation- 
cation interaction. For k ca — 1 and thus T cat — T sta t 
one recovers the standard choice of parameters for the 
CIM. In contrast to previous works we wish to break with 
charge neutrality and rather interpret H sta t as a general 
static disordered energy landscape. Thus T sta t can be 
varied in analogy to o> thereby modifying the strength 
of the disorder. The total Hamiltonian for the CIM is, 



— 1 cat ' — 7} h 1 stat ' ~rp — • (, iU J 



D. Computational Settings 

The number of lattice sites in one dimension was set 
to 20 for all data throughout this work. For simulating 
bulk properties we apply periodic boundary conditions 
and the minimum image convention ]l4| . As shown in 
in Fig. |3| the smallness of the finite size effects justifies 
our choice of 20 lattice sites per dimension. All presented 
data were calculated with a cation concentration of 0.03. 
The number of anions for generating the static energy 
landscape in the CIM was kept identical to the number 
of cations. The contributions of the Coulomb interaction 
were calculated via Ewald summation. The number of 
different starting configurations ranges from 5 to 10. This 
rather small number made it possible to simulate a large 
range of parameters. For model systems with c = 0.03 
and e 2 = 1, which were used throughout this paper, one 
has V c = 0.501. 



The diffusion constants are taken from (r 2 (£))/i-data. 
Fig. ^ shows some examples. Each curve begins with a 
short-time plateau changes into a dispersive regime from 
t\ to t% and forms a long-time plateau beyond £2 . In the 
dispersive regime between t\ and ti the curve follows a 
Jonscher type ]T5| power law 

(r 2 (t))/(6t)~t fc '- 1 (11) 

The following relations are used to determine the diffu- 
sion constants, 

6D(t) = <r 2 (t))/i, 
6D DC =lim^ 00 ,(r 2 (i))/i, 
6D AC =lim t ^ {r 2 {t))/t. 

The high frequency diffusion constant Dac was also cal- 
culated from the average hopping rate a 2 {w) = QDac- 
For details on the general behavior of the REM and the 
CIM, see || H @. 

In the REM the presence of cation-cation interaction 
leads to a significant decrease in diffusion and an increase 
in dispersion Jl3[ . A quantitative analysis of the behav- 
ior can be found in |l7| . For a REM without long range 
Coulomb interaction, i.e. T cat = 0, the activation energy 
E a can be calculated as the difference between the Fermi 
energy and the critical percolation energy |?], [llj. For a 
concentration of 0.03 the Fermi energy is — 1.88<r e and 
the percolation energy is — 0.49er c . Hence the activation 
energy is — 0.49er e + 1.88cr £ = 1.39cr e . This reasoning is 
only valid for low temperatures ay 3> 1. For high tem- 
peratures oy <C 1 it can be shown |l7j that the activation 
energy is roughly equal to o e j \pH = 0.56cr e . Here we ana- 
lyze the REM including long range Coulomb interaction 
for which this analytical treatment is no longer possi- 
ble. Simulated data for various T cat and ay are shown in 
Fig. |^(left). The data for T cat — correspond to vanish- 
ing cation-cation interaction and the regression for the 
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Figure 6: Dependency (CIM) of the low frequency diffusion 
constant Ddc on T ca t- The maximum (T ca t,max) shifts to 
larger T ca t for higher values of F a t a t- The errors are smaller 
than the symbols. 
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Figure 7: Dependency (CIM) of the high frequency diffusion 
constant Dac on T C at- The ) shifts to 

larger T cat for higher values of T s tat- The errors are smaller 
than the symbols. 



low temperature regime leads to an activation energy of 
1.41er e which is in good agreement with the theoretical 
value resulting from percolation arguments, see above. 
One can easily see that the low temperature part of all 
curves can be approximated fairly well by straight lines. 
The first data point for a r = is omitted for all T cat . 
Interestingly, the change of slope at high temperatures, 
i.e. small predicted for T cat = 0, seems to hold 

for large cation-cation interaction, too. The data shown 
in Fig. H(left) can be fitted with the activation energy 
dependent on ay, T cat and a cross term T cat ■ oy, 



In [6 ■ Ddc] = — {a\a r + a 2 T cat + a 3 <T r T cat ) 



(12) 



The slope of the fits in Fig. [|(left) can be expressed 
according to eqn. ( |l2| ) as ai + a3r cat . Fig. fright) il- 
lustrates this dependency. The three coefficients turn 
out to be ax = 1.41, a 2 = 0.063 and a 3 = 0.015. The 
non-vanishing value of a 3 is of particular interest since 
it contradicts the conclusions of Maass et al. |l7). They 
showed that for c=0.01 the low temperature behavior is 
activated, i.e. ln[6Z)] oc ^ when varying the temperature. 
Since ay oc ^ and r cat oc ^ this statement is identical to 
a 3 = 0. We have included a curve in Fig. ^(left) which 
represents a variation of the temperature only, this curve 
shows a proportionality to A/T + B/T 2 (A,B > 0). In 
contrast to this result experimental data for low tem- 
peratures show no curvature |17|] or opposite curvature 
p8[ . This temperature dependence is also very differ- 
ent to that observed in Coulomb glass for which, if at 
all, a description with B < would be appropriate for a 
limited temperature regime in order to recover the Efros- 
Shklovskii law. 

Furthermore we also analyzed the dependency of Dac 
on T sta t and r cat . Fig. ^ shows data for Dac i n anal- 
ogy to those in Fig. [|(left). Exponential fits arc not as 
accurate as for Ddc- Though the data cannot be accu- 
rately fitted as straight lines an approximated slope is 
decreasing with increasing T cat . Thus Dac has a differ- 
ent behavior as compared to Ddc- 
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Figure 8: Correlation between static potential and mean value 
of the dynamic potential for the CIM (left; r cat = 10; T at at = 
10) and the REM (right; F ca t = 40, ov = 4). Each site is 
represented by a single diamond symbol (20 3 for each graph). 
In the CIM a strong correlation is observed, the slope is close 
to one. In the REM no significant correlation is observed. 



The analysis for the CIM was performed in a similar 
way. The first surprising feature is the increase in mo- 
bility of the cations for Ddc (Fig- and Dac (Fig. 
^) when increasing T cat from with fixed T sta t- This 
is the completely opposite behavior as compared to the 
REM. Increasing T cat further leads to an increase in dif- 
fusion until a maximum is reached, from this maximum 
the diffusion decreases and the slope becomes constant in 
a logarithmic plot. Obviously, a simple functional form 
as eqn. ( p"2| ) can not be found for the CIM. A differ- 
ence between Ddc and Dac is that T ca t >max for Ddc is 
approximately equal to T stat whereas T cat ^ max for D A c 
increases more rapidly with r stat . 

The different behavior of the REM as compared to the 
CIM can be rationalized by a closer inspection of the dif- 
ferences in the potential surfaces. Comparing Fig. [j] and 
Fig. ^ and taking the definitions of the models into ac- 
count it's evident that the static potential of a lattice site 
in the CIM-surface is spatially correlated, whereas in the 
REM no correlation among adjacent sites exists. This 
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Figure 9: Matrices of successful and unsuccessful trial jumps 
in the CIM (top, Y cat = Y ata t = 10) and the REM (bottom, 
a r —4, r cat =40) resolved for static and dynamic jump ener- 
gies, z is a counter of how many jumps happen per energy 
interval during a simulation run. The choice of parameters 
has no influence on the qualitative picture. 



correlation in the CIM has a direct consequence: Fig. || 
illustrates a correlation of the static potential of a lat- 
tice site with the mean value of the dynamic potential. 
For a single starting configuration the static energy of 
each site is constant during simulation, being denoted as 
Estatic/kBT in Fig. ||. The mean value of the dynamic 
energy Ed yrm mic a t some site is generated by averaging 
over all energies which a particle at this site has due to 
Heat during a simulation. The observed correlation in 
Fig. ||(left) indicates a correlation of H cat and H an ; high 
cation-cation energy corresponds to low static energy and 
vice versa. It could be shown that the three regions ob- 
served in the CIM correspond to three different types of 
lattice sites. There are sites with no, one or two adjacent 
anions, i.e. anions with a distance of -s/3 • §. No such 
correlation is observed for the REM. 

To elucidate the impact of this correlation on the ion 
dynamics all occurring trial jumps, accepted or not by the 
Metropolis algorithm, were recorded during a MC simu- 
lation with respect to their static (AE stat i C ) and dynamic 

namic 

) jump energy differences (Fig. ^). The sum 
of AE static and AE dynamic is the total energy difference 
for a jump. One apparent effect of the correlation for 
the CIM is quite obvious: a high energy contribution 
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Figure 10: Radial density function g(r) for the CIM(left) and 
the REM(right). For the CIM the parameters are chosen 
to show the temperature dependency and for the REM the 
dependency of g(r) on T cat . The errors are within the size of 
the symbols. 



from the static energy is accompanied by a low energy 
contribution from the dynamic energy and vice versa. 
Therefore to first approximation the total energy does 
not change during a jump giving rise to faster dynamics. 
The REM does not show this behavior. For the REM 
both contributions are independent and both impede the 
cation dynamics. A consequence is a different jump en- 
ergy distribution in the CIM as compared to the REM. 
To further explore the reasons for this behavior the radial 
density function g(r) of both models is illustrated in Fig. 
|lO| . Obviously, the cations in the REM have a structure, 
they prefer a certain distance to each other which returns 
periodically in the graph. This structure is very similar 
to that found in simple liquid systems. The CIM does 
not show such a structure, instead it shows an unusual 
high concentration of cations in very short distances and 
beyond r ~ 3 the bulk density is already reached. The 
found structures give possible explanations for the ob- 
served behaviors: In the REM each cation is surrounded 
by a cage of other cations, and this cage slows down the 
movement p9| because most moves would increase H cat . 
In the CIM the cations prefer to populate low energy 
positions around the anions, which results in an effective 
screening of the anion charge. This screening flattens the 
overall energy landscape and thus gives rise to increased 
mobility. Increasing T cat reduces the clustering due to 
cation-cation repulsion but also supports screening due 
to higher cation charge. These two effects compete and 
may account for the observed behavior, see Fig. ^| and 
Fig. 0. 

We are now able to reason the parameter dependency 
of eqn. ([l2|) for the REM. The first term arises from the 
static disorder of the energy landscape. The second term 
is introduced by a cage effect as mentioned above. The 
cross term on the other hand has no single microscopic 
origin for Doc and Dac- For Due it may arise from 
relaxation of the system to a small perturbation (i.e. a 
cation jump). By a jump process a cation may have left a 
well-adjusted environment and is surrounded by an ener- 
getically unfavorable environment. Apart from jumping 
back the total system may also adjust to this new situ- 
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ation as already formulated in the concept of mismatch 
and relaxation by Funke (2[ ||] . The situation at the new 
position improves energetically after some time due to the 
subsequent relaxation of the adjacent particles, thereby 
reducing the probability of the backjump with time. In 
the presence of disorder this neighbor relaxation is, of 
course, much slower since also the neighbors experience 
the effect of the static disorder. Therefore it is more 
likely for the central particle to jump back, giving rise to 
a decrease of D DC due to the simultaneous effect of static 
disorder and the cation-cation interaction. For Dac an 
increase in T cat reduces the dependency on o~ r and there- 
fore the above argumentation fails. For no cation inter- 
action and in the limit of zero temperature the cations 
occupy only the N lowest sites in static energy. By intro- 
ducing cation interaction the emerging liquid structure 
forces the particles to occupy also some sites which are 
higher in static energy (but lower in total energy), this 
would lead to the observed reduction in a r dependence. 

IV. SUMMARY 

The REM and the CIM display major qualitative dif- 
ferences in their dependences of Dac and Ddc on the 
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